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Abstract 

A  Rouse  model  for  polymer  chains  is  incorporated  into  the  linear  continuous  stick- 
slip  molecular-based  tube  reptation  ideas  of  Doi-Edwards  and  Johnson-Stacer.  This 
treats  the  physically  constrained  (PC)  molecular  stretches  as  internal  strain  variables 
for  the  overall  PC/chemically  cross-linked  (CC)  system.  It  yields  an  explicit  system  of 
stress-strain  equations  for  the  system  permitting  simple  calculations  of  complex  stress- 
strain  relations.  The  model  that  is  developed  here  treats  PC  molecule  as  entrapped 
within  a  constraining  tube,  which  is  comprised  of  both  CC  and  PC  molecules.  The 
model  is  compared  with  experimental  data  sets  from  the  literature. 


1  Introduction 

One  of  the  most  widely  used  empirical  models  for  viscoelasticity  in  materials  is  the  Boltzmann 
convolution  law  [15,  19,  20,  39];  for  a  nice  summary  and  further  references  see  Chapter  2 
of  [20].  In  recent  literature  [5,  8,  9],  models  for  hysteretic  damping  in  elastomers  entail  a 
phenomenological  Boltzmann-type  constitutive  law  of  the  form 

/■*  d 

a{t)  =  ge{e{t))  +  Coeit)  +  Y {t  -  s)  —  {e{s) ,  e{s))  ds  (1) 

J-oo  ds 

where  Y  is  the  convolution  memory  kernel,  and  ge  and  g^  are  nonlinear  functions  accounting 
for  the  elastic  and  viscoelastic  responses  of  the  elastomers,  respectively.  Previous  efforts 
summarized  in  [6]  have  shown,  through  comparison  with  experimental  data,  that  the  best 
ht  to  hlled  elastomer  data  occurs  when  g^  and  gy  are  cubic,  along  with  F  as  a  distribution  of 
exponentials.  Banks,  et  al.,  [7,  2]  subsequently  developed  nonlinear  models  based  on  stick-slip 
“molecular”  ideas  of  Johnson  and  Stacer  [26]  and  Doi  and  Edwards  [10]  which  resulted  in  a 
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form  for  g^,  Qv  and  Y  that  matched  the  empirical  hndings  reported  in  [8,  9,  6].  These  models 
allow  for  multiple  relaxation  times  present  in  polymer  strands  of  composite  materials  within 
a  virtual  compartmental  model  of  entangled  chemically  cross-linked/physically  constrained 
system  of  long  chain  “molecules”.  While  accounting  for  multiple  relaxation  parameters,  the 
models  do  not  include  physically  or  chemically  based  parameters  in  the  polymer  strands. 

In  the  present  paper  a  new  model  is  developed  which  combines  the  virtual  stick-slip  con¬ 
tinuum  “molecular-based”  ideas  of  Johnson  and  Stacer  [26]  with  the  Rouse  molecular-bead 
ideas  as  described  in  Doi  and  Edwards  [10].  This  new  model,  in  which  polymer  chains  are 
treated  as  strings  of  interconnected  beads,  permits  the  incorporation  of  many  important 
physical  parameters  (such  as  temperature,  segment  bond  length,  internal  friction,  and  seg¬ 
ment  density)  in  the  overall  hysteretic  constitutive  relationship.  Our  goal  here  is  to  present 
development  of  this  model  based  upon  physical  considerations  at  the  molecular  level;  its 
form  is  similar  to  that  developed  in  [6,  7]  and  does  have  the  general  form  (1)  of  Boltzmann 
type,  even  though  the  kernel  is  not  of  convolution  type,. 


2  Description  of  the  Rouse  Model 

For  our  summary  of  the  Rouse  model  for  free  polymer  strands  and  subsequent  stress  cal¬ 
culations,  we  follow  for  the  most  part  the  development  in  Doi  and  Edwards  [10],  modi¬ 
fying  somewhat  the  random  noise  assumptions  and  the  particular  series  used  in  order  to 
insure  convergence.  We  hrst  assume  a  material  composed  of  free  polymer  strands  with  each 
strand  consisting  of  a  hnite  set  of  beads  connected  in  a  string  with  elastic  springs.  Let 
{Ri,  R2,  ■■■,  Rn)  be  the  position  vectors  relative  to  a  hxed  coordinate  system  of  the  beads 
comprising  an  interconnected  chain  as  depicted  schematically  in  Figure  1.  Moreover,  let  the 
equation  of  motion  of  the  beads  be  described  by  the  Langevin  equation  [10] 

(9  /  dll  \  1  (9 

+  /„(«)  )  +  (2) 

\  UiLYn  /  ^  OriyYi 


where  fm{t)  is  a  random  force  term,  ks  is  Boltzmann’s  constant,  and  T  is  the  temperature. 
The  mobility  tensor  and  the  interaction  potential  are  chosen  to  be 


c 


^  n=2 

respectively,  with 

SksT 

where  b  is  the  segment  bond  length  at  equilibrium  and  C  is  the  friction  constant  of  the 
polymer  sample. 
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Figure  1:  Representation  of  vectors  for  a  bead-spring  polymer  molecule. 


If  we  use  the  parameters  defined  above  for  the  mobility  tensor,  Hnmi  and  for  the  inter¬ 
action  potential,  U,  then  equation  (2),  for  the  cases  when  n  =  2,3, N  —  1,  can  be  written 
as 

d  R 

=  -k{2Rn  -  Rn+l  -  Rn-l)  +  U  (4) 

at 

For  the  special  cases  of  the  ends  of  the  polymer,  i.e.,  the  cases  when  n  =  1  and  n  =  N,  we 
see  that  (respectively) 


C^  =  -k{RN-RN-l)+fN. 
at 

We  dehne  (R)  to  be  the  conhgurational  average  of  the  beads,  i.e.. 


J  Aip{R;t)dR, 


where  il){R-,t)  is  the  conhgurational  distribution  [10]  of  the  beads.  The  conhgurational  dis¬ 
tribution  is  a  probability  distribution  which  represents  the  probability  that  particles  exist  at 
the  points  (Ri, ...,  R^q)  at  the  given  time  t. 

The  term  /„  =  (/,[,  /^,  is  a  randomly  distributed  force  which  accounts  for  the  Brow¬ 
nian  motion  of  the  beads.  We  assume  therefore  that  the  random  force  /„  is  distribnted 
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according  to  a  Gaussian  distribution  which  is  determined  by  the  following  moments 

=  0, 

{/“(«)/.?(*'))  =  -  «'),  (5) 

where  5ij  =  1  ii  i  =  j  and  0  otherwise,  and  5{t)  is  the  usual  Dirac  function. 

To  obtain  the  equation  of  a  typical  polymer  strand  from  the  hnite  bead  strings,  we 
conceptually  take  a  continuum  limit,  replacing  the  system  (4)  with  an  equation  on  the 
interval  0  <  n  <  N,  where  now  N  is  the  length  of  the  strand.  In  the  limit  we  obtain  a 
partial  differential  equation  in  n  for  the  position  i?(f,  n)  of  particles  along  the  strand  given 
by 


dR(t,  n) 

^  Wt  ^ 

n)  I 
dn 


+  /(*.«).  forO<„<JV 


drR 
dR{t,  n) 
dn 


\n=N 


=  0. 


(6) 

(7) 


This  is  obtained  under  the  assumption  that  Rq  =  Ri  and  Rn+i  =  Rn  and  results  from 
viewing  the  hrst  term  in  the  right  side  of  (4)  as  a  difference  quotient  for  the  second  derivative. 
The  generalized  random  force  f{t,n)  is  now  assumed  to  satisfy 

=  0, 

(/"(f,n)/^(f',m))  =  2(kBTS{n  -  m)Sai3S{t  -  t').  (8) 


A  standard  method  for  analyzing  systems  such  as  (6)-(7)  is  via  Fourier  series  with 
“modal”  coordinates  corresponding  to  the  time  dependent  Fourier  coefficients.  A  system 
of  decoupled  ordinary  differential  equations  is  obtained  through  separation  of  variables  tech¬ 
niques.  That  is,  we  assume  that  R(t,  n)  can  be  expanded  as 

/—  oo 

—  cos(^) 

P=1 


in  terms  of  the  normalized  Fourier  elements  <Pp{n) 
the  random  noise  has  the  form 


We  further  assume  that 


f{t,n)  =  £^pWa/^cos(^) 

p=i 

=  — cos(^),  (9) 

p=i 


where  the  Wp(t)  are  Gaussian  processes  satisfying 


{WM  =  0, 
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(10) 


{u-;(()nf  (('))  =  -  *')■ 

The  coefficients  {/ip}  are  chosen  as 


2  _  12(kBTN _ ^  ^ 

lip  2  2  ’  ’  ’  1 

^  TT^P^ 


(11) 


so  that  all  the  inhnite  series  in  our  subsequent  discussions  below  converge  and  so  that  the 
relationship  in  (8)  is  satished.  We  then  hud  that  the  modal  coordinates  Xp{t)  are  given  by 


r-N 


fpnn 


i?(t,  n)dn,  p  =  1,  2, 


and  satisfy 

d  ^  ^  ^ 

(12) 

C~^Xp  kpXp  +  Qp, 

where 

(13) 

^p=  NW  P  ’  P=l>2,... 

{9p{t))  =  0 

{9p{t)9q{t'))  =  SpqSa06{t  -  t')pl. 


3  Stress  Calculations 

Once  the  modal  coordinates  Xp  have  been  found  for  the  Rouse  model  for  free  polymer 
strands,  it  is  possible  to  use  them  to  determine  a  formula  to  approximate  the  stress  tensor 
for  a  viscoelastic  polymer  undergoing  deformations.  We  use  the  equation  for  the  polymer 
dependent  stress  as  given  by  equation  (7.81)  in  [10],  which  is 

OO 

=  V  E  ( W) 

p=l 

where  c  is  the  segment  “density”  (and  thus  ^  represents  the  number  per  unit  volume  of 
polymer  strands  in  the  solid). 

We  dehne  i?(— 0,n)  to  be  the  position  of  the  polymer  segment  before  deformation,  and 
R(+0,n)  to  be  the  position  of  the  segment  immediately  after  deformation.  Thus  under  the 
affine  deformation  assumption  (e.g.,  see  p.  241,  [10])  which  is  a  linearization  approximation 
(  see  p.  112,  [10]) 

R(+0, n)  =  E(0)  ■  R(-0, n),  (15) 

or  in  terms  of  the  modal  coordinates, 

Xp(+0)  =  E(0)-Xp(-0),  (16) 
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where  the  tensor  E(0)  =  {Eq,^(0)}  is  the  usual  conhguration  gradient  at  time  t  =  0. 

The  matrix  E  is  sometimes  called  the  deformation  gradient  (in  a  misnomer)  but  the  actual 
deformation  gradient  is  D  =  E  —  I  (see  for  example  [4,  5,  8,  29,  31]).  We  recall  that  E  can 
be  used  to  dehne  the  Green-St.  Venant  strain  S  =  ^(E^E  —  I)  =  |(D^D  +  D  +  D"^)  as  well 
as  the  left  Cauchy-Green  strain  =  EE^  which  is  the  same  as  the  Finger  strain  dehned 
below.  The  equation  (12)  for  X^(t)  can  be  solved  to  obtain 


Xfit)  =  A';(0)e-4‘  +  1  y  9<‘(s)e5<-‘>rfs, 


(17) 


where 


tr  = 


STT^ksT 

is  the  Rouse  relaxation  time  (p.  96,  196  of  [10]).  If  we  multiply  on  both  sides  of  (17), 


we  have 


E^^X^it)  =  +  f  '^ds. 

s  Jo 


In  a  similar  manner  we  see  that 


{E^,x';(t)Ef,x;(t))  =  E^^Ef.{x;{t)x;{t)}.  (is) 

Noting  from  (9)  that  gp(t)  =  fipWpit)  and  using  (10),  we  hud  that  the  equation  for  the 
autocorrelation  function  is  given  by 

(X^(t)X^(t)}  =  (A“(0+)A'»(0+)>e-=4‘  +  ^1  _  .  (19) 

We  assume  that  the  system  is  at  equilibrium  before  the  initial  deformation  at  time  t  =  0.  We 
further  assume  that  as  t  approaches  inhnity,  the  system  will  return  to  its  initial  conhguration 
(the  equilibrium  state).  That  is. 


Mm  X;(t)  =  A,“(0-), 

t^OO  ^  ‘ 

which  implies 

(a7(o-)a;(o-)>  = 

If  the  linearization  approximation  (16)  is  used,  we  hud 


where 


(A“(0+)A'’(0+))  =  B„g(E(0))X 

3 

R,^(E(0))  =  J2Eap{O)E0M  =  (E(O)E^(O)),^ 

M=1 


(20) 


(21) 
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is  the  Finger  strain  (see  p.  242,  [10]).  We  note  that  the  Finger  strain  is  related  to  the 
Green-St.  Venant  strains  through  Bap{E?^)  =  2£oti3  +  <^a/3- 

We  may  now  substitute  equation  (20)  into  Equation  (19),  and  we  find  that 


5a/3(E(0)));;^e  + 


2CK 


2Ckp 


JafS 


-2l^t 

e 


(22) 


From  this  expression  we  may  note  that  without  the  fip  terms  as  defined  in  (11),  the  resulting 
series  for  the  stress  tensor  given  by  (14)  does  not  converge!  Also  note  equation  (22)  holds 
for  t  >  0  small  (see  [10]). 


4  Connection  with  “Stick-Slip” 

Equation  (22)  coupled  with  (14)  describes  the  contribution  of  each  node  of  a  free  long  chain 
polymer  molecule  to  the  molecule’s  overall  stress.  The  goal  of  this  work,  however,  is  not  to 
simply  reproduce  a  stress-strain  law  based  on  a  molecule  in  free  space  (i.e.,  based  on  the 
Langevin  equation),  but  rather  to  describe  the  stress  of  a  system  composed  of  physically 
constrained  (PC)  molecules  entangled  with  chemically  cross-linked  (CC)  molecules  and  ex¬ 
periencing  the  stick-slip  mechanisms.  More  precisely,  we  will  view  conceptually  the  material 
undergoing  deformation  as  composed  of  two  virtual  compartments  as  depicted  in  Figure  2. 
One  compartment  will  consist  of  a  constraining  tube  which  is  a  macroscopic  compartment 


Eiih  apped  Portion  of  PC  molecrrle 


Escaped  Poitioii  of  PC  molecule 


Figure  2:  PC  molecule  entrapped  by  the  surrounding  constraining  tube. 

containing  both  CC  and  PC  molecules.  The  other  compartment  will  be  microscopic  in  na¬ 
ture  and  consist  of  those  PC  molecules  aligned  with  the  direction  of  the  deformation.  These 
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molecules  will  at  first  “stick”  to  the  constraining  tube  and  be  carried  along  with  its  motion, 
but  will  very  quickly  “slip”  and  begin  to  “relax”  back  to  a  conhguration  of  lower  strain 
energy.  We  wish  to  compute  the  contributions  of  both  “compartments”  to  the  overall  stress 
of  a  polymer  material  undergoing  deformations. 

To  accomplish  this  goal  we  must  consider  the  contribution  from  the  constraining  tube 
composed  of  both  non-aligned  physically  constrained  molecules  and  chemically  cross-linked 
molecules,  and  that  of  PC  molecules  aligned  in  the  direction  of  the  deformation  that  are 
initially  entangled  with  molecules  of  the  tube.  These  aligned  molecules  will  in  time  escape 
entanglement  and  become  “free”  molecules  and  will  thus  contribute  to  the  overall  stress  in 
two  distinct  phases:  when  entrapped  and  after  “leaking”  free.  Therefore,  there  are  three 
contributions  to  the  stress  of  the  system  the  PC  chain  in  entanglement,  the  portion 

of  the  PC  chain  that  has  escaped  entanglement,  and  the  contribution  due  to  the  constraining 
tube.  The  constraining  tube  will  be  treated  as  elastic  while  we  use  the  Rouse  formulation  to 
treat  the  aligned  PC  molecules.  We  will  denote  the  stress  of  the  portion  of  the  polymer  chain 
that  is  constrained  by  the  surrounding  molecules  as  and  the  stress  of  the  portion  of 

the  polymer  chain  that  has  leaked  out  of  the  constraint  tube  as  The  total  stress 

contribution  of  the  entangled  PC  molecules  will  be  denoted  as 

We  will  denote  the  stress  of  the  constraining  tube,  assumed  to  be  elastic,  by  Thus, 

the  total  polymer  dependent  stress  will  be  formulated  as 

d?(i)  =  d”’(«)  +  d?”’(i).  (23) 

We  will  use  the  Rouse  model  (14)  in  conjunction  with  a  step-strain  process  (similar  to  the 
stick-slip  molecular  formulation  of  Johnson  and  Stacer  [26])  to  arrive  at  an  appropriate  form 
for  u^[g^(t).  This  will  result  in  a  hysteretic  term  as  in  a  Boltzmann-type  stress-strain  law. 
To  calculate  the  contribution  of  the  entangled  portion  of  the  molecule  to  the  stress  we  will 
subject  the  molecule  to  a  series  of  instantaneous  step-strains  at  times  0  =  fo,  •  •  • ,  tn  with 
At  =  ti  —  ti-i  very  small  and  investigate  the  behavior  of  the  component  <  Xp(t)X^(t)  > 
after  each  step-strain,  where  the  PC  molecules  remaining  in  the  tube  are  momentarily  free 
and  thus  subject  to  the  Rouse  dynamics.  We  make,  of  course,  the  additional  assumption 
that  during  the  successive  deformations  some  of  the  entrapped  molecules  will  ’’leak”  and 
escape  entrapment  (this  ’’leaked”  portion  of  the  molecules  will  then  be  considered  free,  so 
we  will  therefore  also  assume  the  Rouse-like  expression  (22)  to  describe  their  motion).  To 
treat  the  entrapped  molecules  we  will  let  the  time  between  each  succeeding  step-strain  go 
to  zero  in  order  to  obtain  a  constitutive  law  that  describes  each  node’s  contribution  to  the 
stress  after  an  instantaneous  step-strain  is  applied  to  the  molecules.  To  arrive  at  that  point. 


first  note  in  (22)  that 

A  (X‘(0)X^(0))  ^  (X“(0+)X^(0+))  -  (X^(O-)X^(O-)) 

(x-(0+)x^(0+))  (a;«(0+)x«(0+)) 

^«/3(E(0))4; 

_  -Ba/3(E(0))  —  Sa0 

-Ba/3(E(0)) 

A(i?.^(E(0))) 

Bafs(E(0))  ^  ^ 

where  Ai?Q,/3(E(0))  =  -Bq,/3(E(0))  —  Sa/3  with  Sa/3  representing  the  Finger  strain  in  the  unde- 
formed  state.  A  simple  manipulation  yields 

(a,"(o+)a;*’(o+)>  =  (A;(o-)A;f(o-)>  +  (b„3(e(o)))  . 

Since  the  PC  molecule  behaves  according  to  Rouse’s  model  momentarily  after  an  instanta¬ 
neous  deformation  (when  it  is  still  considered  free),  (19)  implies 

|(ATWA'fW>  =  ^  . 

or,  equivalently, 

and  we  determine  on  a  short  time  interval  0  =  to  A  ^ 

2  2 

(A»(()A''(«)>  =  + 

=  (A7(0-)A^^(0-)>  +  Ce=^'.  (25) 

According  to  (25)  just  after  the  step  strain  at  t  =  O’*" 

(a;(o+)a;>+)>  =  (A;“(o-)A^^(o-)>  +  Ce"^"* 

=  (A7(0-)A^^(0-))+C,  (26) 

which,  according  to  (24)  implies 


C 


(A7(0+)Af>(0+)> 

(B««(E(0))) 


A(B„ME(0))) 
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and 


(X‘(0-)X^(0-)} 

(A7(0+)X«(0+)) 

(Bo^(E(0))) 


A(B„4E(0)))e 


-2p^ 

-^R 


(t-O) 


(27) 


ioT  0  =  to  <  t  <  ti-  The  above  procedure  will  serve  as  a  basis  for  imitating  Johnson  and 
Stacer’s  step-strain  procedure.  In  order  to  do  so,  we  now  need  to  make  an  assumption  similar 
to  (24)  at  ti.  That  is,  we  would  expect 


or 


A  {X‘(0)Xl{h)) 

{x^(tt)x^(tt)) 


(x-(tt)Xl{tt)) 

AB„f,(E(f,)) 


(A7((pA«((+)> 


\^P  v^i  // 

ABa/3(E(ti)) 


(A';(«r)A^^(i+)) , 


where 

Ai?Q,/3(E(ti))  =  Ba/3(E(ti  -h  At))  —  Ba/3(E(t^  )). 
Thus,  it  is  necessary  to  understand  what  the  quantity 


(28) 


(29) 


(■V(y)y''(y)>  -  (A';(«-)Ay(i-)> 

in  general  represents.  To  investigate  this  quantity  we  compute 


(A';(«,  +  A«)A^(«,  +  At))  -  {X‘(t-)Xl(t-)) 

for  At  >  0  and  then  let  At  tend  to  zero  from  above.  This  procedure  leads  to  the  conclusion 
{X^iU  +  At)x^^(t,  +  At))  -  {X^{t-)Xl{t-))  ^  [AB^f,{E{U))/At]  At 


where 

We  remark  that 


ABai3{E(ti))  —  Bai3(E{ti  J-  At))  —  Bais(E(ti)). 


lim 

At^0+ 


AB^fs{E{ti)) 

At 


5S„,(E(t)) 


t=U 
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may  either  exist  in  the  ordinary  sense  or  may  provide  a  jump  aX  t  =  ti.  Arguments  for  this 
approximation  and  further  details  on  the  quantity 


(X‘(u  +  Xt)X^(U  +  A«)>  -  {X‘(t-)Xl(t-)) 


are  found  in  the  appendix. 

Continuing  with  our  arguments  based  on  the  Johnson  and  Stacer  step-strain  procedure, 
we  recall  from  (29)  and  (30)  that 


{x;(tt)xii(tt)) 


{Xi'(0-)x^(0-))  + 

A  (A7(y)Ay(y)) 


A(B„^(E(i,)))eA<‘ 


-ti) 


Since  Rouse’s  model  requires  (X"(f)X^(f))  (X"(0— )X^(0— ))  exponentially,  we  hnd 


(A';(o-)yyo-)> 


(  A  {x;(tt)xntt)) 


A{Bap(E{ti)))e  -R 


e 


for  ti  <t  <t2-  If  this  process  is  repeated  indehnitely,  we  hnd  that 


{X‘(t)X^^(t)) 


(A7(0-)A<'(0-)> 

A  AAdJAil) 

A  (B„«(E(i.))) 


A(B„j(E(ii)))eA<‘  '■> 


(31) 


for  t  >  tk- 

As  mentioned  above,  we  assume  that  some  portion  of  the  entrapped  molecule  can  escape 
and  behave  as  a  free  molecule.  To  account  for  these  dynamics,  we  dehne  ■yit)  to  be  the 
fraction  of  the  molecule  that  is  still  entrapped  at  time  t  (so  that  7(0)  =  1).  Recalling  that 
N  is  the  length  of  the  molecule,  we  dehne  Ne{t)  =  y{t)N  as  the  length  of  the  molecule  still 
entrapped  at  time  t.  Thus  for  the  entrapped  portion  we  have  that 

2(kp 


Returning  to  (22)  we  hnd  if  we  let  N  =  Ne  then  the  entrapped  molecule  contributes 


252  iV3 


(XyijX^t)}  = 


=^(t-to) 


2b^N^ 


T{t)  1  —  e 


(32) 
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(33) 


to  the  stress  for  0  =  to  <  ^  <  ^i,  which  leads  to  the  approximation 

{X^{4)X^{4))  2b^N^ 


TT'^P'^ 


-7  (fo 


-Bo/3(E(to)) 

Since  the  relaxation  of  the  PC  molecules  obeys  Rouses’s  model  for  a  very  short  time  after 
the  instantaneous  step-strain,  it  follows  that  on  ti  <  t  <  tj+i,  (32)  holds  with  to  replaced  by 
ti-  This  leads  immediately  to 

(X“(t+)Xp^(t+)>  2b^N^ 


TT'^P^ 


■7  (ti) 


for  ti  <t  <  tj+i-  Doi  and  Edwards  (p.  196,  [10]  or  [21])  calculate 


7(f)  = 


e 


where 


Td  = 


p  odd 

CN%‘^ 


Tr'^ksTa'^ 


is  the  disengagement  time.  If  (34)  is  substituted  into  (31)  we  hnd 

{x^{t)x^{t))  «  (x;(o-)x^^(o-)> 


^  2b^N^  , 

E 


i=0 


7^  (t,) 


At 


At, 


for  t  >  tk-  Taking  the  limit  as  At  =  t  —  ti  goes  to  zero,  we  obtain 

(V»(«)A'»(«)>  =  (X“(0-)X^(0-)) 

2b‘^N^  /■*  n,  .  d  ,,,  =^(t-s)  , 

7^(s)— (Ra/3(E(s)))e  "fl  ’ds. 


ttV  Jo 


Therefore,  the  contribution  to  the  stress  of  the  constrained  molecule  is  given  by 


(34) 


Aw  = 


N 


P=1 


P=1 


/i. 


2Ck 


^  5 


c  ^  ,  2b^N 

+ 


vrV 


E 

p=i 

oo 

z 

P=1 


GcksT 

7j-2p2 

QcksT 

jj-2p2 


dai3+  I  7^('S)^  (5„^(E(s)))e  -R  "Us 


j  E(t,s)^(R„^(E(s)))ds 
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where  V (t,  s)  =  7^(s)F (t  —  s)  with 


^ ,  ,  QckBT  ^  1 

p=l  "  p=l 


E' 


=  >  CpC 


(35) 


for 


Cp  — 


QabsT 

P^2p2 


and  np  =  ^. 


We  note  that  this  stress  term  can  be  written  in  terms  of  the  left  Cauchy-Green  strain  as 


d 


p=i 


^2p2 


ds 


(36) 


Observe  that  the  form  of  o'd)(t)  is  the  similar  to  that  of  the  Boltzmann-type  stress-strain 
law  (1)  except  that  the  kernel  Y  is  no  longer  in  simple  convolution  form  as  in  (35). 

The  contribution  to  the  stress  due  to  the  portion  of  each  polymer  chain  that  has  leaked 
out  of  the  constraint  tube  is  given  by  a  modihcation  of  (22)  similar  to  the  one  used  in  (32) 
above 


-2p^ 


^2Ck 

p=i  ' 


^  ^  Kn/u  '^0/3(1  e 


^2CA:p 


E 

p=i 


GchsT 

7j-2p2 


-2p^ 


(l-7(t))"5„;3(E(0))e 


^  QckBT  j  X  /  1  -r- 

WWW  (1  “  ^(^))  (  1  -  e 


p=i 


^2p2 


=^t 


Finally,  for  the  contribution  of  the  elastic  constraining  tube  to  the  overall  stress 

in  (23)  we  choose  the  Finger  strain  Hookean  form 


W’(i)  =  p''S„p(E(i)) 


(37) 


where  is  a  generalized  Young’s  modulus  of  elasticity. 


5  Uniaxial  Deformation 

In  this  section  we  consider  uniaxial  deformations  and  examine  the  equation  for  the  macro¬ 
scopic  stress  [10],  which  is  of  the  form 

W/3(t)  =  +  P6ap.  (38) 
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Here  the  term  represents  the  contribution  from  the  polymer  molecules  (the  polymer 
dependent  stress,  as  dehned  in  equation  (23))  and  P  is  the  hydrostatic  pressure. 

We  proceed  by  assuming  that  we  are  applying  a  tensile  deformation,  i.e.,  a  deformation 
strictly  in  one  of  the  three  principle  directions  (specihcally,  we  will  consider  a  stretch  in  the 
z  =  X3  direction).  To  determine  the  stress  for  such  a  deformation,  hrst  the  appropriate 
Finger  strain  is  required.  If  we  consider  a  unit  cube,  and  apply  a  small  deformation  in  the 
2;  direction,  then  it  attains  a  length  of 

A  =  1  +  e, 

where  the  strain,  e,  is  the  ratio  of  the  change  in  length,  AL,  to  the  original  length,  L,  or 
in  other  words,  e  =  ^.  If  the  material  is  assumed  incompressible,  the  volume  must  be 
maintained.  Thus  the  sides  in  the  x  =  Xi  and  y  =  X2  direction  must  both  necessarily  be  of 
length 

We  choose  a  random  point  within  the  solid  denoted  by  .R  =  The  point’s 

change  in  position  after  deformation  from  R  to  the  new  location,  denoted  by  R,  can  be 
described,  to  hrst  order,  by  the  equations  (p.  241,  [10]) 

pi  _  ^  pi  p2  _  ^  p2  p3  _  \  p3 

it  —  — 1=^  s  -rt  —  — ,  it  —  Ait  . 

These  equations  give  a  conhguration  gradient  E  of  the  form 

0  0] 

0  0  A. 

which  provides  a  Finger  strain  of  the  form 

"A  0  O' 

B(E)=  0  i  0  . 

_0  0  A2_ 

We  dehne  the  tensile  stress  in  the  principle  direction  using  the  macroscopic  stress 
given  by  equation  (38) 

y  =  AP)  I  p 

where  P  is  a  hnite  hydrostatic  pressure  term  [11,  31].  For  our  case  of  deformation  in  the  2; 
direction,  we  consider  the  equation 

E,  =  yfi  +  p, 

in  which  we  must  determine  the  hydrostatic  pressure  term  P.  This  is  done  by  noting  that 
since  the  deformation  is  uniaxial  in  the  2:-direction,  no  force  acts  in  the  x-  or  y-directions. 
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Thus,  the  stress  in  the  x-  and  y-directions  vanishes,  i.e.,  =  0,  and  from  the  equation 

Sx  =  Uxx  +  P  for  the  tensile  stress  in  the  x  direction,  it  is  seen  that  P  =  —aW- 

Thus,  if  a  tensile  deformation  is  performed  on  the  elastomer  along  the  axis  then  the 
stress  is  given  by 


S.(t)  = 


= 


'{t)  +  P 


- 


=  a 


= 


it) 

it)  -  (^xJit)  +  (^fKt)  -  (^xxit)  + 

t 


[elas) 

zz 


it) 


a. 


(elas) 


it) 


QckRT  ^ 


p=i 


TT^ 


1  r  3/  N 
^  /  7“(s)e  'V 
P  Jo 


2X(s)X\s)  + 


A2(.) 


ds 


(1  -7(^))' 

p2 


X\0)  - 


m 


=jpt 


+  X\t)- 


m 


(39) 


The  term  jp  (A^  —  ,  where  jp  is  the  Young’s  modulus  of  elasticity,  accounts  for  the  stress 

contribution  of  the  elastic  constraining  tube.  We  observe  that  this  corresponds  to  the  Cauchy 
or  true  stress  for  an  incompressible  neo-Hookean  material  undergoing  uniaxial  elongation. 
This  can  be  derived  in  a  pseudo-phenomenological  approach  [4,  5,  8]  using  the  Mooney  strain 
energy  function  (SEF)  in  the  context  of  a  nonlinear  elasticity  approach  [29,  31,  34,  38,  39]. 


6  Parameter  Estimation  and  Simulation  Results 

6.1  Articular  cartilage  results 

In  this  subsection,  we  report  on  calculations  performed  with  the  stress-strain  relationship 
(39)  using  parameters  determined  from  a  set  of  data  from  experiments  on  articular  cartilage 
(a  material  of  signihcant  scientihc  interest  which  is  widely  viewed  as  a  viscoelastic  material- 
see  [21]  and  the  references  therein).  This  will  provide  a  first  test  of  our  stress  models 
in  reproducing  results  from  other  models  and  physical  experiments.  First  the  appropriate 
parameters  used  in  the  model  are  estimated  in  inverse  problems  incorporating  the  data.  Then 
stress  calculations  which  are  based  on  the  results  of  experimental  work  will  be  presented. 
The  stress-strain  relation  will  be  evaluated  by  calculating  the  stress  and  comparing  it  to 
experiments  for  various  input  strain  functions.  Finally,  simulations  are  conducted  to  repeat 
the  results  of  Johnson  and  Stacer’s  paper  [26]  on  which  the  model  is  partially  based. 

6.1.1  Estimating  parameters  and  corresponding  stress-strain  simulations 

The  experiments  conducted  by  Huang,  et.ah,  [24]  involved  applying  a  tensile  strain  (defor¬ 
mation)  to  a  sample  of  articular  cartilage  and  then  measuring  the  stress  within  the  cartilage. 
Two  such  experiments  were  conducted  in  which  two  different  input  strains  were  used.  These 
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strains  were  ramp  strains  starting  at  zero  and  increasing  at  a  constant  rate  until  a  cessation 
time  (ts).  The  material  was  then  held  at  a  hxed  strain  emax  until  the  experiment  terminates 
at  time  tj  =  2000  seconds.  For  both  experiments  e^nax  was  taken  to  be  0.05,  while  the  first 
had  a  cessation  time  of  t]  =  .126  seconds  and  the  second  had  a  cessation  time  of  tl  =  400 
seconds.  Thus  the  equation  for  the  strain  functions  is  given  by 


eiit) 


t  <  ti 


^max 


t>f. 


(40) 


for  i  =  1,2.  The  chosen  parameters  used  in  the  strain  function  for  these  experiments  are 
presented  in  Table  1  while  the  graphs  of  the  strain  functions  are  given  in  Figure  3. 


Table  1:  Parameters  used  for  ramp  input  strain  function  for  cartilage  stress  experiments. 


Parameter 

Abbreviation 

Value 

Cessation  Time  1 

tl 

.126  seconds 

Cessation  Time  2 

tl 

400  seconds 

Maximum  Strain 

^max 

0.05 

Strain  versus  Time 


Figure  3:  Input  strains  used  for  stress  calculations  in  cartilage  stress  simulation. 

It  is  assumed  that  the  experiments  were  conducted  at  room  temperature  which  is  taken 
to  be  T  =  300°  K.  For  the  length  of  the  polymer  chain  we  chose  N  =  1000.  A  sensitivity 
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analysis  of  the  model  (39)  with  respect  to  parameters  such  as  the  step  length  a,  bond  length 
b,  segment  density  c,  the  frictional  constant  (  and  the  constraining  tube’s  elastic  constant 
jj,  =  jjX  was  carried  out  to  determine  possible  correlations.  Based  on  these  investigations  the 
model  was  re-parameterized  by  dehning  a  =  ^  (the  parameter  characterizing  r^),  b  = 

Y 

(the  parameter  characterizing  r^j),  and  /i  =  ^  (a  normalized  Young’s  constant).  (We  remark 
that  in  view  of  (43)  and  (44)  below,  we  could  have  included  and  in  the  scaled  variables 
b  and  a,  respectively,  but  the  estimation  results  would  have  ultimately  produced  the  same 
hts  to  data  in  the  efforts  with  experimental  data  reported  on  below.)  To  obtain  values  for 
the  parameters  a,  b,  c,  and  jl,  parameter  estimation  methods  using  the  experimental  data 
with  the  model  were  employed  . 

To  perform  the  parameter  estimation,  data  was  extracted  from  the  graphs  presented  in 
Figure  7  of  [24].  These  graphs  depict  the  stress  on  articular  cartilage  for  applied  ramp  strains 
as  described  above.  (Graphs  of  the  extracted  data  are  presented  as  solid  lines  in  Figure  4.) 
The  data  was  extracted  using  the  MatLab  tool  Grabit,  written  by  Jiro  Doke  [12].  Two  sets 
of  data  were  obtained,  one  from  each  of  the  experiments  performed,  and  these  are  referred 
to  as  and  where  is  the  stress  value  for  the  strain  function  eiitj)  as 

described  by  equation  (40))  at  time  tj. 

The  data  is  then  used  in  a  weighted  least-squares  cost  function  to  determine  the  optimal 
values  for  the  desired  parameters  in  a  vector 


r  0,  1  r  a  ] 


L  94  J  L  A 1 


The  weighted  least  squares  (WLS)  function  is  given  by 

2  .  N  2  r 100  2' 

(42) 

with  which  we  employed  a  Nelder-Mead  method  {fminsearch  in  MatLab)  to  determine  the 
optimal  value  of  the  parameter  vector  6.  In  this  case  all  of  the  parameters  we  seek  are 
uniquely  dehned.  The  response  function  is  dehned  by  letting 

Aj(s)  =  2Ai(s)A'(s)  -h 

and  then  dehning  the  function 
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The  summation  limit  M  was  chosen  to  be  10  since  it  was  found  that  the  total  sum  would 
change  by  less  than  1%  if  M  were  increased.  In  computing  the  integrals  in  (43),  an  approx¬ 
imation  was  made  using  the  MatLab  quad  routine,  which  utilizes  a  Simpson’s  quadrature 
method  to  approximate  the  integral.  We  also  used 


(44) 


where  the  limit  of  summation  was  chosen  to  be  M'  =  21  (for  a  larger  M'  value,  the  increase  in 
the  sum  is  not  sufficient  to  justify  the  increased  computation  time).  The  term  Aj(t)  =  H-ej(t) 
is  dehned  for  ej(t),  as  given  above.  The  derivative  of  \i{t)  was  taken  for  i  =  1,2,  as 


ei(t) 


^max 

0 


t<tl 

t>t\ 


The  cost  function  was  calculated  using  equation  (43)  in  the  minimization  algorithms.  These 
calculations  require  an  initial  guess,  denoted  by  6*o,  for  the  value  of  6.  Since  little  can  be  found 
in  the  literature  for  these  parameters  in  the  case  of  cartilage,  in  this  section  we  obtained 
initial  values  Oq  by  simulating  with  the  model  with  numerous  parameter  values  over  a  wide 
range  and  comparing  the  corresponding  graphs  visibly  with  the  data.  A  physically-based 
method  for  determining  initial  estimates  in  the  case  where  one  knows  rough  parameter  ranges 
for  a  material  is  described  for  the  case  of  polyisoprene  data  in  the  next  section. 

In  addition  to  calculating  the  optimal  values  of  the  parameters  in  9  we  will  determine  the 
standard  errors  for  a,  b  and  fi.  The  process  to  calculate  the  standard  errors  depends  upon 
the  form  of  the  cost  functional.  If  the  ordinary  least  squares  (OLS)  cost  functional  is  used 
(as  it  will  later  in  (45))  then  the  standard  error  is  approximated  as 


SEu{9^)  =  v'Cfcfc(r) 


The  term  Ckk  is  the  diagonal  element  of  the  M  x  M  covariance  matrix 


where 


c{n  =  cr^[x^{hx{h 
dm 


-1 


Xjk{9)  = 


d9k 


is  the  (j,  k)  element  of  the  sensitivity  matrix  y  G  ,  9^  G  3?^  is  the  parameter  estimate 
obtained  in  the  optimization  process,  and 


d2  = 


n  —  M 


E  /.r)-» 


i=i 
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with  fjiO"^)  =  T,z{tj',9^)  denoting  the  model’s  value  at  time  tj  and  parameter  estimate 
6^.  More  details  regarding  large  sample  size  approximation  statistics  can  be  found  in  the 
standard  nonlinear  regression  approximation  theory  ([13,  22,  25],  and  Chapter  12  of  [35]). 
For  a  brief  summary  also  see  Section  3  of  [3]. 

If,  on  the  other  hand,  a  weighted  least  squares  cost  functional  is  used  (as  in  (42))  a  matrix 
that  accounts  for  the  individual  weights  must  be  included  in  the  formation  of  the  covariance 
matrix 

r  .  .  .  n -1 


C{9^) 


=  a 


X^{9^)W{9^)X{9^) 


The  weights  used  in  (42)  produce  the  weighting  matrix 


W  \9)  =  diag  (max{|/^}, . . . ,  max{|/^},  max{|/^}, . . . ,  max{|/^}) 


We  note  that  the  segment  density  c  acts  as  a  scaling  parameter  for  the  model  and  hence 
requires  a  single  data  point  to  set  its  value.  Therefore,  analysis  reveals  that  the  model  is 
relatively  insensitive  to  changes  in  c,  making  a  standard  error  calculation  involving  multiple 
observations  irrelevant. 

For  the  case  when  9o  is  taken  to  be 

7  X  10-6  ■ 

9  X  10“^ 

10-3 
10^ 


the  optimal  value  returned  by  the  MatLab  program  is 


3.4592  X  10-6  ±  1.1828  x  10-^6 
4.0786  X  10-3  ±  3.3634  x  lO-^^ 
5.7831  X  10-^ 

1.5587  X  10^  ±  1.1383  x  10-^3 


These  optimal  values,  along  with  the  other  important  values  associated  with  the  calcu¬ 
lation  of  the  stress  function,  are  presented  in  Table  2. 


Table  2:  Fixed  and  optimal  parameter  values  used  for  cartilage  stress  simulations 


Parameter 

Abbreviation 

Value 

Temperature 

T 

300°  K 

Segments  /  Chain 

N 

1000 

Boltzmann’s  Constant 

ks 

1.3806505  X  10-23  J/K 

Disengagement  constant 

^opt 

3.4592  X  10-6  A^kg^ 

Rouse  constant 

bopt 

4.0786  X  10-3  A^kg/s 

Segments  /Volume 

^opt 

5.7831  X  10-^  1/A3 

Young’s  constant 

f^opt 

1.5587  X  102  MPa  A3 
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Once  we  have  determined  a  set  of  parameters  which  provide  an  optimal  £t  to  the  data 
obtained  from  the  paper  by  Huang,  et  ah,  we  now  can  perform  calculations  which  will 
simulate  the  experiments. 

The  resulting  simulations  are  presented  in  Figure  4  in  a  comparison  with  the  data  ob¬ 
tained  from  [24], 


Stress  versus  Time  Stress  versus  Time 


Time,  t  (s)  Time,  t  (s) 

Cartilage  Experiment  1  Cartilage  Experiment  2 

Figure  4:  Stress  calculation  with  optimal  fit  parameters  versus  experimental  data  for  cartilage 
with  ramp  strain  inputs. 


Instead  of  using  both  sets  of  data  to  obtain  9opt  as  we  did  in  minimizing  (42),  we  could 
also  obtain  a  set  of  optimal  parameters  for  each  experiment  separately.  That  is,  we  could 
use 

100  2 

=  ,  (45) 

i=i 

to  obtain  separate  optimal  parameters  9lp^,i  =  1,2,  for  each  experiment.  When  the  experi¬ 
ments  are  considered  separately,  the  corresponding  solutions  with  might  better  approx¬ 
imate  the  data  from  experiment  i  than  those  with  9opt  obtained  using  (42). 

For  the  data  from  the  first  experiment  we  set 

1  X  10“^ 

8  X  10“^ 

2  X  10-3 
4 

and  after  optimization  obtained 

3.6712  X  10-3  ±  1.3978  x  lO-^^ 

7.3624  X  10-^  ±  5.8073  x  lO-^^ 

3.3313  X  10-3 

4.4014  X  10^  ±  3.5449  x  lO-^^ 
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With  the  second  data  set,  we  used  6*q  =  6*g  and  obtained 

■  3.6955  X  10“^  ±  1.6159  x  lO'^^  ' 

^  _  4.0242  X  10-3  ±  8.4251  x  lO'^^ 

“  4.6266  X  10“^ 

1.9261  X  lO^  ±  2.7588  x  lO'^^ 

A  comparison  of  corresponding  stresses  with  6lp^  with  the  data  for  each  experiment  is  pre¬ 
sented  in  Figure  5.  It  is  obvious  from  the  second  experiment  that  optimizing  with  its  data 
set  separately  produces  parameters  that  yield  a  model  that  more  closely  agrees  with  the 
data;  the  results  from  the  hrst  data  set  reflect  the  model’s  improved  ability  to  achieve  the 
data’s  peak,  but  in  doing  so  a  portion  of  the  tail  of  the  data  is  missed. 


Stress  versus  Time  Stress  versus  Time 


Time,  t  (s)  Time,  t  (s) 

Cartilage  Experiment  1  Cartilage  Experiment  2 

Figure  5:  Stress  calculations  and  experimental  data  for  cartilage  with  ramp  strains  using 
separate  cost  functions  (45). 

6.1.2  Stress  versus  strain  model  simulations 

Having  estimated  parameters  for  the  stress-strain  model  (39),  one  can  then  use  this  model  in 
simulations  with  various  input  strains  to  investigate  the  possible  presence  of  features  such  as 
hysteresis.  For  example,  when  a  stress-strain  relation  appears  to  possess  a  simple  one-to-one 
graph  in  response  to  periodically  oscillatory  strain  inputs,  this  indicates  that  there  is  little 
strain  rate  dependence  in  the  system,  i.e.,  no  hysteresis  is  present.  However,  one  expects 
that  the  graph  of  the  stress-strain  relation  will  appear  as  loops  in  response  to  such  inputs 
when  the  system  contains  hysteresis. 

In  a  series  of  simulations,  various  strain  functions  were  input  to  the  stress  function  of 
(39),  and  the  stress-strain  relations  for  each  were  graphed.  For  each  of  these  simulations 
in  this  subsection,  the  parameter  values  used  were  those  given  in  Table  2;  the  strain  input 
functions  vs.  time,  the  resulting  stress  vs.  time  and  the  stress  vs.  strain  relation  are  plotted 
below.  For  each  simulation,  the  input  strain  function  is  taken  on  the  interval  from  t  =  0 
seconds  to  tf  =  200  seconds. 
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For  the  first  simulation  a  sinusoidal  input  strain  function  given  by 


1  1 Q  y| 

ei(t)  =  emax{-7:  cos(^t)  +  1)  (46) 

Z  tf 

is  used  to  test  the  response  of  the  system  to  cyclical  input.  This  particular  sinusoid  was 
designed  so  that  it  is  never  negative  and  it  ranges  between  0  and  emax  =  0.05.  It  can  be  seen 
from  the  results  presented  in  Figure  6  that  the  stress  exhibits  hysteresis. 


Strain  vs.  Time 


Stress  vs.  time 


Stress  vs.  Strain 


Figure  6:  Graph  of  input  strain  ei(t)  =  emaxi—^  cos{^t)  +  1),  model  stress  vs.  time,  and 
model  stress  vs.  strain  relations  for  ei  and 
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In  the  next  simulation  we  used  sinusoidal  strain  function  with  increasing  amplitude  given 


by 

/  \  /  1  /  IOtt  n  ,  n  t  / , 

=  emax{-^COs{-—t)  +  1)  —  .  (47) 

^  tf  if 

This  simulation  was  performed  in  an  attempt  to  see  if  the  stress  of  the  system  would  vary 
in  a  way  other  than  linearly  with  increasingly  greater  strain  cycles.  The  results  are  depicted 
in  Figure  7. 


strain  versus  Time  Stress  versus  Time 


Strain  vs.  Time  Stress  vs.  time 


Stress  versus  Strain 


Stress  VS.  Strain 

Figure  7:  Graph  of  input  strain  e2(t)  =  emaa;(— |  cos(^f)  +  l)^,  model  stress  vs.  time,  and 
model  stress  vs.  strain  relations  for  62  and 

For  a  third  simulation,  a  simple  bell  curve  for  the  strain  input  is  employed.  This  particular 
strain  function  is  chosen  because  of  its  simplicity.  The  function  used  to  describe  this  strain 
is  dehned  by 

r3W  =  ^(-cos(^«)  +  l).  (48) 

which  is  a  single  period  of  a  cosine  function.  Results  for  this  input  strain  function  are 
presented  in  Figure  8,  in  which  it  may  be  seen  that  the  stress-strain  relation  is  a  simple  loop. 
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strain  versus  Time 


Stress  versus  Time 


Strain  vs.  Time  Stress  vs.  time 


Stress  versus  Strain 


Stress  VS.  Strain 

Figure  8:  Graph  of  input  strain  e^it)  =  cos{j^t)  +  1),  model  stress  vs.  time,  and 

model  stress  vs.  strain  relations  for  €3  and  S^. 
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6.2  Results  for  Polyisoprene 

We  next  investigated  use  of  the  model  (39),  or  equivalently  (43),  with  experimental  data 
for  polyisoprene.  For  the  model  calculations  of  the  stress  for  polyisoprene,  it  is,  of  course, 
necessary  to  determine  the  particular  parameters  associated  with  that  polymer  (specifically 
the  parameters  a,  b,  c,  N,  and  /i  must  be  estimated). 

We  will  choose  the  temperature  to  be  at  298°K  which  is  roughly  room  temperature. 
There  have  been  many  experiments  conducted  for  polyisoprene  at  this  temperature,  and 
thus  the  amount  of  data  from  which  we  may  derive  hrst  estimates  of  some  of  the  parameters, 
especially  the  density  parameter  c  and  the  friction  parameter  (  (a  factor  in  both  a  and  b) 
which  are  affected  by  temperature,  is  substantial. 

Using  knowledge  of  the  monomer  chemical  structure  (see  Figure  9)  of  polyisoprene,  we 
can  hrst  calculate  a  parameter  Xp  referred  to  as  the  degree  of  polymerization  which  is  simply 
the  number  of  monomers  per  single  molecule.  The  degree  of  polymerization  can  be  calculated 
using  the  equation 

M 

with  molecular  weight  M  and  the  monomer  weight  Mq  (also  referred  to  as  the  mass  of  the 
repeat  unit).  We  can  calculate  the  monomer  weight  using  values  observed  in  the  monomer 
structure.  We  have  that  per  monomer  there  are  5  carbon  atoms  and  8  hydrogen  atoms. 

TCH, 

^C=C 

H3C  H 


Figure  9:  Chemical  structure  of  polyisoprene  polymer. 

Thus,  if  we  obtain  the  atomic  weight  for  carbon  (which  is  12.01)  and  for  hydrogen  (which  is 
1.08)  from  any  standard  periodic  table,  we  have  that 

Mo  =  5(12.01)  +  8(1.008) 

=  68.114. 

For  the  molecular  weight  M,  we  use  the  value  presented  by  Stille  (Table  10.2,  [36])  which  is 
M  =  350,  000.  We  are  then  able  to  calculate  the  degree  of  polymerization  Xp  to  be 

350,000 
“  68.114 
«  5138, 

where  we  have  rounded  to  the  nearest  whole  monomer.  The  bond  length  b  is  determined 
by  segmentation  of  the  polymer  where  the  polymer  chain  is  segmented  into  N  segments  of 
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length  b.  These  are  related  to  the  contour  length  L  hj  L  =  Nb.  As  argued  in  [23],  b  can 
be  approximated  by  6  ~  8.415  A  while  L  =  Xpbm  where  bm  is  the  average  monomer  length 
given  approximately  [37]  by  4.602  A  so  that  L  ^  23647.25  A.  It  then  follows  from  N  =  L/b 
that  N  ^  2810. 

We  also  need  to  determine  an  approximate  value  for  the  segment  density  c  of  the  sample. 
This  is  calculated  using  the  formula  (see  equation  (2)  in  [27]), 

_  pNa 

Mo  ’ 

where  p  is  the  polymer  density,  Na  is  Avogadro’s  number,  and  Mq  is  the  monomer  molecular 
weight.  The  value  of  Mq  was  given  above  as  Mq  =  68.114.  For  the  polymer  density  we  use 
the  value  given  by  Abdel-Goad,  et  ah,  [1],  as  p  =  0.90  x  10“^^  g/A^  at  T  =  298  K  (values 
range  from  0.9  x  10“^"^  to  1  x  10“^"^  g/A^,  [14], [16], [17], [18],  but  most  commonly  0.9  x  10“^^ 
g/A^).  The  segment  density  is  then 

_  rNa 

Mo 

_  0.90  X  10-24(6.02  X  1023) 

“  68.114 

~  0.007954  segments/A3. 

We  consider  next  the  friction  constant  C,.  In  the  Ferry  text  ([15]),  the  friction  coefficient 
C  is  shown  to  be  equal  to  the  product  of  the  number  Xp  of  monomers  per  molecule  (the 
degree  of  polymerization),  and  the  monomeric  friction  coefficient  C,m  (the  friction  provided 
by  a  single  monomer).  In  Table  12-III  on  page  258  of  [15],  the  monomer  friction  coefficient 
for  polyisoprene  (listed  as  unvulcanized  Hevea  rubber)  at  room  temperature,  T  =  298°  K,  is 
given  in  log  form  to  be 

log(Cm)  =  —6.74  dynes  s/cm. 

When  converted  to  our  units  we  obtain 

Cm  =  1.18265  X  10“®  kg/s. 

Then  when  we  apply  the  formula  for  C  we  hnd  that 

C  =  5138Cm 
~  0.006076  kg/s 

is  the  friction  present  in  a  single  strand  of  the  polymer. 

All  that  remains  is  the  step-length  a.  Doi  and  Edwards  [10]  dehne  a  in  their  equation 
(6.4)  by 

_ 

—  ) 
a 
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where  L  is  the  contour  length  (which  was  calculated  earlier  to  be  approximately  23647.25 
A)  and  N  is  the  number  of  segments  of  length  b.  (We  note  that  Nb"^  is  the  mean  square 
end-to-end  distance  (r^)o  of  the  chain  [10]).  We  thus  hnd  a  =  8.415  A.  We  can  then  use 
these  values  of  a,  b  and  (  to  provide  initial  estimates  for  a  and  b.  We  collect  these  and  other 
pertinent  parameters  used  for  the  polyisoprene  estimation  procedures  in  Table  3. 


Table  3:  Fixed  parameters  along  with  initial  values  used  for  polyisoprene  stress  versus  strain 
estimation. _ 


Parameter 

Abbreviation 

Value 

Temperature 

T 

298°K 

Segments  /  Chain 

N 

2810 

Boltzmann’s  Constant 

ks 

1.3806505  X  10-23  J/K 

Disengagement  constant 

ho 

4.3026  X  10-3  A^kgA 

Rouse  constant 

bo 

4.3026  X  10-3  A^kg/s 

Segments  /Volume 

Co 

7.9540  X  10-3  1/A3 

Young’s  constant 

ho 

1.2572  X  10-3  MPa  A3 

6.2.1  Polyisoprene  parameter  estimation  resnlts 


Having  determined  a  set  of  approximate  parameters  related  to  natural  rubber  (polyisoprene), 
we  then  used  these  as  an  initial  guess  in  an  OLS  estimation  procedure  for  the  “best”  param¬ 
eters  in  (39).  The  data  used  in  the  OLS  was  provided  by  a  graph  from  the  text  by  Riande, 
et  ah,  [33].  We  used  this  data  to  define  an  input  strain  function  in  the  form  of  a  simple  hat 
function 


e{t) 


f-t  t<% 

t>%  ' 


(49) 


as  graphed  in  Figure  10,  with  a  maximum  value  of  emax  =  2  and  tf  =  500  seconds.  We  used 


Strain  versus  Time 


Figure  10:  Strain  function  from  which  the  stress  for  the  polyisoprene  simulations  are  carried 
out. 
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this  hat  function  as  input  for  stick-slip  model  (39)  and  computed  the  corresponding  stress 
to  use  as  the  model  in  an  OLS  cost  functional 


100 


C(«)  =  |E,(«_,i  »)-!,> 

i=i 


(50) 


where  the  stress  data  points  obtained  from  the  graph  from  [33]  are  referred  to  as 
with  each  representing  the  stress  for  e{tj)  for  tj  the  uniformly  spaced  time  points  on  the 
interval  from  t  =  0  to  t  =  500  seconds.  This  stress  function  T,z(tj-,9)  =  of  (43) 

and  the  stress  data  was  used  to  estimate  the  parameters  in  the  model.  In  (43)  we  used 
Xi(t)  =  X(t)  =  1  +  e(t)  for  e  given  in  (49),  and  we  again  choose  the  limits  of  summation  to  be 
M  =  10  and  M'  =  21.  We  carried  out  estimation  of  the  parameters  using  the  initial  values 
for  a,  b,  c  and  /x  given  in  Table  3.  The  optimal  values  (with  standard  errors)  found  are 

2.9746  X  10-1  ±  4.2663  x  lO'^^  ' 

5.7572  X  10-1  ±  1.6687  x  lO-i® 

5.2495  X  10-^ 

2.2017  X  10-2  ±  1.0104  X  lO-i® 


A  graph  of  the  corresponding  model  stress-strain  curve  using  the  optimal  parameter  values 
is  compared  to  the  experimental  data  in  Figure  11.  We  note  that  while  the  basic  scale 


Stress  versus  Strain  Results  for  Polyisoprene  Rubber  -  optimizaed  parameters 


Stick-Slip/Rouse  Hybrid  Model  and  Comparative  Data  from  [33] 

Figure  11:  Model  stress  versus  strain  simulation  using  optimal  parameters  compared  to  data 
for  polyisoprene. 

and  trends  of  the  two  graphs  are  qualitatively  similar,  it  is  interesting  to  note  that  there  is 
little  to  no  hysteresis  exhibited  by  pure  natural  rubber,  both  in  our  simulations  and  in  the 
comparative  data.  Moreover,  there  are  nonlinear  aspects  of  the  data  that  clearly  are  not 
captured  by  the  model.  We  recall  the  linearization  assumption  of  (15)  which  might  suggest 
difficulties  for  the  model  when  used  with  nonlinear  materials. 


6.3  Polyisoprene  with  carbon  black  reinforcement 

Natural  rubber  (polyisoprene),  as  seen  in  the  graph  of  Figure  11,  exhibits  mild  nonlinear 
behavior  but  very  little  hysteresis.  However,  most  rubber-based  products  contain  rubber 
composites  (or  filled  rubber)  and  it  has  been  known  for  a  long  time  that  various  properties 
of  the  composite  are  affected  by  incorporating  substances  (a  common  practice  for  industrial 
products)  such  as  carbon  black  or  colloidal  carbon  into  the  polyisoprene  (i.e.,  by  vulcanizing 
it).  In  the  paper  by  Parkinson  [32],  it  is  noted  that  carbon  black  particles  are  typically 
spherical  and  range  in  diameter  from  roughly  50  to  5000  A,  although  other  sources  have  it  as 
being  between  10  to  10000  A.  When  introduced  into  the  raw  polymer,  the  polymer  strands 
attach  to  these  spheres,  restricting  the  flow  of  the  polymers. 

One  side  effect  of  the  addition  of  carbon  black  to  polyisoprene  is  the  production  of 
significant  hysteresis  in  the  stress-strain  relation  of  the  material.  In  addition,  other,  more 
desirable,  effects  including  an  increase  in  the  stiffness  of  the  material,  resistance  to  absorption 
of  other  fluids,  abrasion  resistance  and  heat  resistance  in  the  composite  may  be  produced. 
These  effects  can  often  be  tailored  to  the  desired  levels  by  varying  the  amount  and  type  of 
carbon  black  introduced  to  the  raw  polymer.  For  example,  carbon  black  reinforced  rubbers 
are  used  in  both  car  tires  and  in  rubber  stoppers.  At  hrst  glance,  these  items  may  appear 
to  be  made  of  different  material,  because  tire  rubber  is  so  much  stiffer  than  a  typical  rubber 
stopper;  however  they  only  differ  in  their  carbon  black  content. 

In  the  text  by  Riande,  et  ah,  [33],  there  is  an  experimental  data  set  for  rubber  reinforced 
with  carbon  black  (graphed  here  as  the  solid  curve  in  Figure  13  below)  undergoing  deforma¬ 
tions.  While  the  concentration  of  the  carbon  black  within  the  sample  used  in  the  experiment 
presented  in  that  graph  is  not  known,  it  is  possible,  if  we  use  the  same  inverse  problem 
methods  as  those  of  subsections  6.1.1,  to  estimate  the  parameters  in  our  model  (43). 

The  data  was  hrst  extracted  from  the  graph  of  [33],  in  a  manner  similar  to  that  used  in 
section  6.1.1.  It  was  assumed  that  the  experiment  takes  place  on  the  interval  0  <  t  <  500 
and  that  the  strain  function  (which  is  not  known  to  us)  was  approximately  piecewise  linear. 
Under  this  assumption  an  approximate  strain  function,  obtained  from  the  data  set,  was 
found  by  determining  the  lines  which  contain  the  maximum  and  minimum  points  of  the 
strain  and  dehning  the  strain  piecewise.  The  strain  function  is  then  of  the  form 

(  (2.07051/165)t  0  <  t  <  170 

eRcsit)  =  <  -(1.9257/165)(t  -  330)  +  0.1448  170  <  t  <  330  .  (51) 

[  (1.98671/170)(t- 330) +  0.1448  t  >  330 

The  graph  of  this  strain  function  is  given  in  Figure  12. 

The  stress  data  points  obtained  from  the  graph  from  [33]  are  referred  to  as 
where  each  is  represents  the  stress  for  eRcsitj)  for  tj  the  uniformly  spaced  time  points 
on  the  interval  [0,500].  This  strain  function  and  the  stress  data  were  used  to  estimate  the 
parameters  in  the  model  (43).  We  assume  that  we  are  working  with  approximately  the  same 
number  of  segments  as  from  Section  6.2.1  and  take  N  =  2000;  we  use  the  same  temperature, 
T  =  298°  K. 
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strain  versus  Time 


Figure  12:  Strain  function  from  which  the  stress  for  the  polyisoprene  with  carbon  black 
reinforcement  simulations  are  carried  out. 


As  we  did  earlier,  we  need  to  find  the  parameters  a,  b,  c,  and  /i.  We  again  dehne  the 
vector  6  of  parameters  to  be 


e  = 


a 

b 

c 

fi 


The  function  A  is  given  by  A(t)  =  1  +  eRCBit),  where  eRcsit)  is  dehned  as  above,  with  a 
corresponding  piecewise  constant  derivative  X'.  The  cost  function  used  in  the  Nelder-Mead 
algorithm  was  the  same  as  that  given  in  (50). 

If  the  initial  guess  is  taken  to  be 


1.0  X  10-3  ■ 
1.0  X  10-3 
1.0  X  10-3  ’ 

10 


then  the  optimal  value  found  is 


1.5556  X  10-^  ±  4.9263  x  lO-^^  ' 

6.5350  X  10-^  ±  9.2709  x  lO-^® 

3.4526  X  10-3 

2.2248  X  10^  ±  1.1351  x  lO-^^ 

The  hxed  and  optimal  parameter  values  used  in  the  model  simulation  for  the  polyiso¬ 
prene  with  carbon  black  reinforcement  are  collected  in  Table  4.  The  model  simulations  are 
compared  to  the  experimental  data  from  [33]  in  Figure  13. 

As  can  be  seen  from  the  graph,  the  model  does  exhibit  hysteresis  in  the  polyisoprene 
simulations  (as  evidenced  by  the  presence  of  loops  in  the  stress-strain  curve).  However, 
the  model  does  not  capture  very  well  the  nonlinearities  in  the  data.  We  recall  that  our 
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Table  4:  Optimal  parameters  used  for  polyisoprene  with  carbon  black  reinforcement  stress 
versus  strain  simulation. _ 


Parameter 

Abbreviation 

Value 

Temperature 

T 

298°  K 

Segments  /  Chain 

N 

2000 

Boltzmann’s  Constant 

ks 

1.3806505  X  10-23  J/K 

Disengagement  constant 

^opt 

1.5556  X  10“^  A^kg^ 

Rouse  constant 

bopt 

6.5350  X  10-^  A^kg/s 

Segments  /Volume 

^opt 

3.4526  X  10“^  1/A3 

Young’s  constant 

f^opt 

2.2248  X  102  MPa  A3 

model  is  based  on  the  linear  Rouse  model  (4)  (which  is  linear  due  to  the  assumptions  on 
the  underlying  potentials)  and  the  linearization  assumption  (15).  Moreover,  both  (24)  and 
(28)  are  assumptions  that  are  linear  in  nature  (as  opposed  to  the  nonlinear  Johnson-Stacer 
type  ratio  assumptions  of  [6,  7]  which  led  to  full  nonlinear  models  for  hysteretic  constitutive 
laws  there).  To  provide  a  better  £t  to  the  carbon-black  hlled  polyisoprene  data,  one  may 
require  nonlinear  assumptions  in  the  model  development  similar  to  those  in  [6,  7],  or  the 
more  appropriate  assumption  in  (2)  of  a  Lennard-Jones  type  higher  order  potential  [28,  30] 
which  permits  repulsive  forces  as  “beads”  in  the  molecular  chain  become  close.  On  a  more 
positive  note,  it  is  clear  that  it  is  possible  to  portray  the  qualitative  trends  in  hysteretic  and 
non-hysteretic  materials  with  the  present  model. 


Stress  versus  Strain  Results  for  RGB 


Stick-Slip/Rouse  Hybrid  Model 

Figure  13:  Stress  versus  strain  calculation  and  comparative  data  for  polyisoprene  with  carbon 
black  reinforcement  simulation. 
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7  Conclusions 

In  this  paper  we  have  derived  a  new  molecular  based  model  to  describe  hysteresis  in  polymers 
and  polymer-like  materials.  The  models,  based  on  a  linear  stick-slip  assumption,  while 
capturing  some  of  the  hysteresis  in  several  types  of  materials,  do  not  capture  especially  well 
the  shapes  of  the  nonlinearities  in  the  hysteresis  loops.  While  these  models  are  a  useful 
initial  effort,  further  development  along  the  lines  of  the  nonlinear  stick-slip  assumptions  of 
[6,  7]  is  needed. 

8  Appendix 

In  this  section  Einstein  notation  (summation  over  repeated  indices)  will  be  used,  so  the 
expression  '^^Eaf^(ti)E/s^(ti)  (our  previous  notation)  is  the  same  as  Ea^{ti)Ei3^{ti)  . 

We  have  claimed  that  the  quantity 

{X;(U  +  +  A())  -  (V“(«r)A«((-)> 

can  be  approximated  by 

-[5«^(E(t))] 

To  verify  this  claim,  we  will  consider  two  cases:  the  configuration  gradient  Ea^{ti)  is  discon¬ 
tinuous  (case  1)  and  the  configuration  gradient  Eafj,(ti)  is  continuous  (case  2). 

Let  Ea^iti)  =  6a/_i  +  and  let  be  the  current  configuration  of  the  polymer. 

Case  1  can  be  resolved  by  observing  that 

+  A()A«((j  +  A())  -  (V“(«r)A«((-)> 

can  be  written,  using  arguments  based  on  those  presented  in  Section  3  and  grouping  higher 
order  terms  (H.O.T.),  as 

=  C{AfSa{tt)  +  A^f,{ti)}+E.O.T. 

=  +  H{At)A^,{ti)][60,  +  H{At)A^,{tt)] 

-[6^^  +  H{0)A^^{ti)][6fs^  +  H{0)Afs^{ti)}  At  +  H.O.T. 

=  {Ea^iti  -\-  At)Ep^(ti  -\-  At)  —  Ea^(ti)Ei3f^(ti)}  At  -\-  H.O.T. 

~  C {Ea^{ti  +  At)Ef^^[ti  +  At)  —  Ea^{ti)Ei^f^[ti)}  At  (52) 

where  C*  is  a  constant  and  the  function  H{x)  is  the  Heaviside  function  H{x)  =  0  for  a;  <  0 
and  H{x)  =  1  for  a;  >  0.  Observe  in  (52)  that 

lim  At)Efj^{ti  -\-  At)  —  Eafj,{ti)Efj^{ti)}  ^  (i?Q^(E(tj  -|-  At))  —  BafjiEiti))) 

At^o+  At  At^o+  At 
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This  implies 


(A';(«+)X|(i+)>  -  (X‘(t-)X^^(t-))  «  c|b„s(E(«))|^_^  Ai. 

For  case  2  we  will  suppose  Ea^(t)  is  differentiable  from  the  right  at  t  =  ti.  Thus 
exists,  Eaf,(tf)  =  limAi^o+  +  At)  and  +  At)  ^  6^^,  +  E'^^(ti)At  (5^^  denotes 

the  configuration  before  the  deformation).  Again,  if  arguments  based  on  those  presented  in 
Section  3  are  used,  then 

(A7((.  +  +  A())  -  (X^(t-)X^(tn} 

can  be  approximated  by 

C  (E'^^iU)  +  E'^JU))  At  +  H.O.T.  =  C  {E'^^{U)E0^{U)  +  E'p^{U)E^,{U))  At  +  H.O.T. 

=  At  +  H.O.T. 

The  above  arguments  justify  the  approximation 

{X‘{u  +  At)X^(U  +  At))  -  (A7((-)A;f'(«-)>  «  Cj^B,,ft(m){_^  At 

in  both  the  continuous  and  discontinuous  cases. 
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